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To efficiently implement many-particle quantum simulations on quantum computers we develop and 
present methods for inverting the Campbell-Baker-Hausdorff lemma to 3rd and 4th order in the com- 
mutator. That is, we reexpress exp {— i (Hi + H2 + . . .) At} as a product of factors exp (— IHiAt), 
exp (—iBi2At), . . . which is accurate to 3rd or 4th order in At. 

PACS numbers: 03.67.Lx 
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Quantum computers have generated much interest re- 
cently, largely due to the result by Shor that they can 
factor integers in an amount of time that grows polyno- 
mially with the size of the integer. This can be compared 
to factorization on a classical computer, where the time 
it takes to factor a number grows exponentially with the 
input size. In addition to Shor's factorization algorithm, 
simulations of quantum systems have also been shown to 
be possible in polynomial time Indeed, this was the 
first area for which it was proposed that quantum com- 
puters could fundamentally be more powerful (i.e. much 
faster) than classical computers || . 

From a theoretical standpoint, a quantum computer is 
a quantum system with a 2"-dimensional Hilbert space. 
Pairs of states in the system are defined to be 'qubits'. 
The canonical example of such a system is a set of n 
spins. Each spin consists of two states, so each spin can 
represent a qubit and the Hilbert space of the system 
is 2™-dimensional. The equivalent of a logical gate on a 
classical computer is an operator acting on a set of qubits 
on a quantum computer. 

This letter focuses on a problem which concerns simu- 
lational issues in quantum computation. A simulation of 
a quantum mechanical system on a quantum computer 
consists of applying an operator exp(-iHt) on a set of 
qubits, where H , the Hamiltonian of the system of inter- 
est, is suitably encoded (and discretized) to act on the 
set of qubits. For many-particle systems H is a sum of 
terms. For instance, the Hubbard model Hamiltonian, 
used in the study of high-T c superconductivity, can be 
written IB] as the sum 



(i) 



where Vb is the strength of the potential, and m a is the 
operator for the number of fermions of spin a at site i. 
In the second (kinetic energy) term, the sum indi- 
cates all neighboring pairs of sites, to is the strength of 



the "hopping" , and Cj, 



are annihilation and creation 



operators, respectively, of a fermion at site i and spin a. 
This model gives an example in which a full simulation on 
a classical computer is impossible due to the exponential 
increase in the size of the Hilbert space of the quantum 
system with the number of lattice sites. 



The canonical quantum computer cannot act on all 
spins at once || . Therefore, it becomes necessary to find 
ways of approximating the evolution operator, which is 
the exponential of a sum of operators (with a Hamilto- 
nian such as that in Eq. (Q)) as a product of operators 
each acting on a subspace of the Hilbert space. To second 
order, for instance, we could use the approximation 
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where the e lH ^ At act on a subspace of the Hilbert space. 

To find higher order approximation methods, we want 
to reexpress exp (^2^=1 An) as a product of individ- 
ual exp(A re )'s. In order to do this, we must invert the 
Campbell-Baker-Hausdorff formula. To 5th order, the 
Campbell-Baker-Hausdorff formula is 
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where 
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As a strategy for finding approximation methods, we 
pick a fundamental ordering of the product of exponen- 
tials with parameters allowing for transposes of the en- 
tire product as well as raising all the exponentials in the 
fundamental unit to the same power. By iterating the 
Campbell-Baker-Hausdorff formula, we can get an ex- 
pression for this fundamental unit in terms of a single 
exponential 
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which defines the J5^ in terms of the A n . Here, p is an 
exponent on a, and a label on the matrices B V N . 
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Now combine a succession i = 1, . . . , I of fundamen- 
tal units with parameters a; and ai. Again iterating 
Campbell-Baker-Hausdorff gives 
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The 73^ are generated from the J5^ by commutation. X 
represents a label pq . . .rs where 



npq...rs = r R p [n« 
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B P ^"' rs is of order p + q- 
can take 
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• r + s. Up to 5th order we 



X e {1; 2; 3, 12; 4, 13, 112; 5, 14, 23, 113, 221, 1112} (8) 

These Bf span the space of the B^'s and their commu- 
tators to 5th order and for TV > 2 they are independent. 
The of are defined in terms of and a.; by Eq. (|^). 
Here again, the X's are labels. After some calculation, 
the Campbell-Baker-Hausdorff formula, Eq. (||), gives the 
equations 
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for p = 1, ... ,5, 
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time evolution during part of the method.^] This follows 
immediately from Eq. (^) with p — 3. It can also be 
proved using Eq. (|9|) with p = 3 and p = 4 that 4th order 
methods must have at least two inverses. 

Our basic method to solve Eqs. (pf|i"l|) is to pick values 
of ai and ai and see if they satisfy the equations. To do 
this we must restrict the number of fundamental units by 
fixing /. We also take the ai's to be ±1 and restrict the 
range of the a^'s. 

We start with Eq. ([)]), since, in this equation, order 
with respect to i does not matter. So, for a given set of 
values, we need to consider only one permutation, not all 
permutations of the values. This greatly reduces the size 
of the search. 

Furthermore, we start by considering p = 1 and 3, 
since it is only the sign of o^Oi that matters in these 
equations. This means we can consider only the sign 
of the combination o^eti, and not the signs of ai and 
ai individually. This reduces the search further. These 
equations are particularly restrictive for the case of few 
inverses. 

After solving the p — 1 and 3 equations, we introduce 
separate signs for the a^s and a^'s and solve the equation 
with p = 2, and p — 4 for the 4th order case. 

Finally, into the restricted set of solutions to Eq. (||) we 
introduce permutations of the a^s and a^s with respect 
to the index i and solve Eqs. (p|~[TT|) . 

We find a larger number of solutions than we can easily 
present. We want to present solutions which are in some 
sense optimal. To do this, we consider the form of the 
operator resulting from a given method 
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for ppq = 112, 113, 221, where of = -of 
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For approximations to exp fj^nLi -^n] > we require all 
of = except for o\ which is the coefficient of B X N = 
Yln=i ^n, and which should be greater than zero. 

An interesting feature of 3rd order methods is that they 
require at least one inverse, i.e. they require backward 



n(^ 



-iajAi At — iajA 2 At — ia,jA N At\ a i 



exp 



JV 



io}J2 An &t + r(-i At)° +1 



n=l 



(13) 



where At <C 1 is a time step, o is the order of the method, 
and 
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where X £ {4, 13, 112} for a 3rd order method and X £ 
{5, 14, 23, 113, 221, 1112} for a 4th order method. 

r is an error which takes values in the vector space 
of the commutators for which we do not have a metric. 
Therefore, we make an ad hoc choice of basis to be dis- 
cussed elsewhere. This allows us to replace r by a single 
real scalar R. The error from the method can then be 
taken to be 



* After this work was completed, we became aware that this 
point had also been noted in H . 
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E = nRAt° +1 
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where n is the number of times we apply the approximate 
method. 

If the physical time we want to simulate is T p , then 



T p = nDAt 



(16) 



where D = a\ is given by the method. 

The computer time it takes for a given simulation can 
be written 



T c = nINtg + nLNt s 



(17) 



where / is the number of fundamental units in the 
method and N is the number of terms in a unit, t g is 
the time it takes to make a gate change, 
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so that LN is the total time the gates are applied for 
in the method. The time an individual gate is applied 
for will be t s — b At, where b is a proportionality con- 
stant dictated by the actual couplings in the quantum 
computer hardware. 

Using Eqs. ( p"5| ) and (|l6|), the computer time can be 
rewritten 
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There are two possible limits to this equation. One is 
that the computer time is dominated by gate switching. 
In this case, we want the factor 
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to be small. The second limit is when the computer time 
is dominated by the time during which the gates are ap- 
plied. Here, the ratio L/D should be small, and to make 
the error small, we want (R/ D)(At)° small. However, if 
At can be made very small (from the hardware point of 
view), then making the error small forces the computer 
time to be dominated by gate switching. If there is a 
limit to At, and it is reached before the computer time is 
gate switching dominated, then the computer time may 
still be dominated by gate application and we want L/D 
and R/D small. We also prefer to have concise methods. 

The 3rd order method that we have selected given the 
above criteria is 
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and the 4th order method is 
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To illustrate our methods, we have applied first, sec- 
ond, third and fourth order methods to the exactly solu- 
ble operator 
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where C = cos (-\/3 At) and S = sin (\/3 At) . 

As a measure of the error, we took the differences Aa x , 
Ady and Aa z between the a x , a y and a z components of 
the exact solution and those of the results of our methods. 
We then calculated the error 
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In Fig. (Q), we plot the logarithm of the error as a 
function of the logarithm of the time that the system was 
evolved for. The first order method results are uppermost 
and higher order results lie underneath each other with 
fourth order results being the lowermost plotted. At = 
0.01 for all methods. 




FIG. 1. A measurement of the accuracy of our results. 
Plotted is the log(error) of (from top to bottom) first, second, 
third and fourth order approximation methods as a function 
of log(time). Note that for the fourth order method, the error 
never grows larger than 10 -3 . 
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Notice that the first order error oscillates once it 
reaches order 1. The rest of the errors remain small 
throughout the simulation, with the fourth order error 
remaining below 10 -3 for the entire evolution. 

The error for all methods goes as nR(At)° +1 , where 
n is the number of times the method has been applied. 



Therefore, logE = log n + log R{At)° +1 



For At = 



0.01, this makes the y-intercept decrease roughly by order 
—2 as the order of the method increases. Since the time 
evolved is proportional to n, the slope of the errors is 1 
for all methods. 

As an example of how useful our approximations can 
be, let us consider a case in which we want to apply an 
approximation method for time T = 1 with total error 
E = 10 -4 . For a first order method, this means that 
we require about 5000 applications of the method. For 
second order, we require about 30 applications. For our 
third order method, we need 2 applications. And for our 
fourth order method, we need less than 1 application of 
the method. This results in a reduction of orders of mag- 
nitude in the computational cost of a given simulation. 
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